import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 正常显示负号

# 简化的共形映射实验
print("实验三：共形映射")

# 简化版映射函数
def map_to_half_plane(z):
    # 将单位圆盘与圆盘|z-1/2|<1/2的交集映射到上半平面
    return 1 / (z - 1)

# 创建交集区域的点
r = np.linspace(0, 1, 15)
theta = np.linspace(0, 2*np.pi, 30)
R, Theta = np.meshgrid(r, theta)
unit_disk_points = R * np.exp(1j * Theta)

intersection_points = []
for i in range(unit_disk_points.shape[0]):
    for j in range(unit_disk_points.shape[1]):
        z = unit_disk_points[i, j]
        if abs(z - 0.5) < 0.5:
            intersection_points.append(z)
intersection_points = np.array(intersection_points)

print(f"交集区域包含 {len(intersection_points)} 个点")

# 应用映射
if len(intersection_points) > 0:
    mapped_points = map_to_half_plane(np.array(intersection_points))
    finite_points = mapped_points[np.isfinite(mapped_points) & (np.abs(mapped_points) < 100)]

    # 可视化
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # 原始区域
    unit_circle = np.exp(1j * np.linspace(0, 2*np.pi, 100))
    ax1.plot(np.real(unit_circle), np.imag(unit_circle), 'b-', linewidth=1, label='|z|=1')
    circle2 = 0.5 + 0.5 * np.exp(1j * np.linspace(0, 2*np.pi, 100))
    ax1.plot(np.real(circle2), np.imag(circle2), 'r-', linewidth=1, label='|z-1/2|=1/2')
    ax1.scatter(np.real(intersection_points), np.imag(intersection_points), s=1, alpha=0.6, color='green', label='交集区域')
    ax1.set_xlim(-0.2, 1.2); ax1.set_ylim(-0.7, 0.7); ax1.set_title('原始区域'); ax1.grid(True); ax1.legend()
    
    # 映射后的区域
    ax2.scatter(np.real(finite_points), np.imag(finite_points), s=1, alpha=0.6, color='purple', label='映射后的区域')
    ax2.set_xlabel('实部'); ax2.set_ylabel('虚部'); ax2.set_title('映射到上半平面'); ax2.grid(True); ax2.legend()
    plt.tight_layout(); plt.show()

# 验证共形映射的保角性质
print("\n验证共形映射的保角性质:")

# 创建网格
x = np.linspace(-1.5, 1.5, 10)
y = np.linspace(-1.5, 1.5, 10)
X, Y = np.meshgrid(x, y)
Z = X + 1j*Y

# 共形映射 w = z^2
def conformal_map(z):
    return z**2

W = conformal_map(Z)

# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# 原始网格
ax1.plot(X, Y, color='blue', linewidth=0.5)
ax1.plot(X.T, Y.T, color='blue', linewidth=0.5)
ax1.set_xlim(-2, 2); ax1.set_ylim(-2, 2); ax1.set_title('原始网格'); ax1.grid(True)

# 映射后的网格
ax2.plot(np.real(W), np.imag(W), color='red', linewidth=0.5)
ax2.plot(np.real(W.T), np.imag(W.T), color='red', linewidth=0.5)
ax2.set_xlim(-3, 3); ax2.set_ylim(-3, 3); ax2.set_title('w=z^2 映射后的网格'); ax2.grid(True)

plt.tight_layout(); plt.show()

# 计算角度保持情况
z0 = 1 + 1j
delta = 0.1
v1_original = delta  # 沿实轴方向
v2_original = 1j * delta  # 沿虚轴方向

w0 = conformal_map(z0)
w1 = conformal_map(z0 + v1_original)
w2 = conformal_map(z0 + v2_original)

v1_mapped = w1 - w0
v2_mapped = w2 - w0

angle_original = np.angle(v2_original / v1_original)
angle_mapped = np.angle(v2_mapped / v1_mapped)

print(f"原始向量间角度: {np.degrees(angle_original):.2f}度")
print(f"映射后向量间角度: {np.degrees(angle_mapped):.2f}度")
print(f"角度差: {abs(np.degrees(angle_original - angle_mapped)):.6f}度")